C     ALGORITHM 587, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL. 8, NO. 3,
C     SEP., 1982, P.323.
      SUBROUTINE LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, LSEI  10
     * MODE, WS, IP)
C
C     DIMENSION W(MDW,N+1),PRGOPT(*),X(N),
C     WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
C     ABOVE, K=MAX(MA+MG,N).
C
C     DABSTRACT
C
C     THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES
C     PROBLEM WITH BOTH EQUALITY AND INEQUALITY CONSTRAINTS, AND, IF THE
C     USER REQUESTS, OBTAINS A COVARIANCE MATRIX OF THE SOLUTION
C     PARAMETERS.
C
C     SUPPOSE THERE ARE GIVEN MATRICES E, A AND G OF RESPECTIVE
C     DIMENSIONS ME BY N, MA BY N AND MG BY N, AND VECTORS F, B AND H OF
C     RESPECTIVE LENGTHS ME, MA AND MG.  THIS SUBROUTINE SOLVES THE
C     LINEARLY CONSTRAINED LEAST SQUARES PROBLEM
C
C                   EX = F, (E ME BY N) (EQUATIONS TO BE EXACTLY
C                                       SATISFIED)
C                   AX = B, (A MA BY N) (EQUATIONS TO BE
C                                       APPROXIMATELY SATISFIED,
C                                       LEAST SQUARES SENSE)
C                   GX.GE.H,(G MG BY N) (INEQUALITY CONSTRAINTS)
C
C     THE INEQUALITIES GX.GE.H MEAN THAT EVERY COMPONENT OF THE PRODUCT
C     GX MUST BE .GE. THE CORRESPONDING COMPONENT OF H.
C
C     IN CASE THE EQUALITY CONSTRAINTS CANNOT BE SATISFIED, A
C     GENERALIZED INVERSE SOLUTION RESIDUAL VECTOR LENGTH IS OBTAINED
C     FOR F-EX. THIS IS THE MINIMAL LENGTH POSSIBLE FOR F-EX.
C
C
C     ANY VALUES ME.GE.0, MA.GE.0, OR MG.GE.0 ARE PERMITTED.  THE
C     RANK OF THE MATRIX E IS ESTIMATED DURING THE COMPUTATION. WE CALL
C     THIS VALUE KRANKE. IT IS AN OUTPUT PARAMETER IN IP(1) DEFINED
C     BELOW. USING A GENERALIZED INVERSE SOLUTION OF EX=F, A REDUCED
C     LEAST SQUARES PROBLEM WITH INEQUALITY CONSTRAINTS IS OBTAINED.
C     THE TOLERANCES USED IN THESE TESTS FOR DETERMINING THE RANK
C     OF E AND THE RANK OF THE REDUCED LEAST SQUARES PROBLEM ARE
C     GIVEN IN SANDIA TECH. REPT. SAND 78-1290. THEY CAN BE
C     MODIFIED BY THE USER IF NEW VALUES ARE PROVIDED IN
C     THE OPTION LIST OF THE ARRAY PRGOPT(*).
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START EDITING AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/DASUM/DASUM/,/DDOT/DDOT/,
C     /DNRM2/DNRM2/,/ DSQRT/ DSQRT/,/ DABS/ DABS/,/DMAX1/DMAX1/,
C     /DCOPY/DCOPY/,/DSCAL/DSCAL/,/DAXPY/DAXPY/,/DSWAP/DSWAP/,/D0/D0/,
C     /, DUMMY/,SNGL(DUMMY)/,/DRELPR/DRELPR/
C
C     WRITTEN BY R. J. HANSON AND K. H. HASKELL.  FOR FURTHER MATH.
C     AND ALGORITHMIC DETAILS SEE SANDIA LABORATORIES TECH. REPTS.
C     SAND 77-0552, (1978), SAND 78-1290, (1979), AND
C     MATH. PROGRAMMING, VOL. 21, (1981), P.98-118.
C
C     THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST..
C     W(MDW,N+1),PRGOPT(*),X(N),WS(2*(ME+N)+K+(MG+2)*(N+7)),IP(MG+2*N+2)
C     WHERE K=MAX(MA+MG,N).  THIS ALLOWS FOR A SOLUTION OF A RANGE OF
C     PROBLEMS IN THE GIVEN WORKING SPACE.  THE DIMENSION OF WS(*)
C     GIVEN IS A NECESSARY OVERESTIMATE.  ONCE A PARTICULAR PROBLEM
C     HAS BEEN RUN, THE OUTPUT PARAMETER IP(3) GIVES THE ACTUAL
C     DIMENSION REQUIRED FOR THAT PROBLEM.
C
C     THE PARAMETERS FOR LSEI( ) ARE
C
C     INPUT..
C
C     W(*,*),MDW,   THE ARRAY W(*,*) IS DOUBLY SUBSCRIPTED WITH
C     ME,MA,MG,N    FIRST DIMENSIONING PARAMETER EQUAL TO MDW.
C                   FOR THIS DISCUSSION LET US CALL M = ME+MA+MG.  THEN
C                   MDW MUST SATISFY MDW.GE.M.  THE CONDITION
C                   MDW.LT.M IS AN ERROR.
C
C                   THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS
C
C                                  (E  F)
C                                  (A  B)
C                                  (G  H)
C
C                   IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1
C                   RESPECTIVELY.
C
C                   THE INTEGERS ME, MA, AND MG ARE THE
C                   RESPECTIVE MATRIX ROW DIMENSIONS
C                   OF E, A AND G. EACH MATRIX HAS N COLUMNS.
C
C     PRGOPT(*)    THIS ARRAY IS THE OPTION VECTOR.
C                  IF THE USER IS SATISFIED WITH THE NOMINAL
C                  SUBPROGRAM FEATURES SET
C
C                  PRGOPT(1)=1 (OR PRGOPT(1)=1.0)
C
C                  OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF
C                  GROUPS OF DATA OF THE FOLLOWING FORM
C
C                  LINK
C                  KEY
C                  DATA SET
C
C                  THE PARAMETERS LINK AND KEY ARE EACH ONE WORD.
C                  THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS.
C                  THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY.
C                  THE VALUE OF LINK POINTS TO THE FIRST
C                  ENTRY OF THE NEXT GROUP OF DATA WITHIN
C                  PRGOPT(*).  THE EXCEPTION IS WHEN THERE ARE
C                  NO MORE OPTIONS TO CHANGE.  IN THAT
C                  CASE LINK=1 AND THE VALUES KEY AND DATA SET
C                  ARE NOT REFERENCED. THE GENERAL LAYOUT OF
C                  PRGOPT(*) IS AS FOLLOWS.
C
C               ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP)
C               .  PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE)
C               .       .
C               .       .
C               .       .
C               ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF
C               .                       NEXT GROUP)
C               .  PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(LINK1+2)=DATA VALUE
C               ...     .
C               .       .
C               .       .
C               ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS.
C                  A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR.
C                  THIS HELPS PREVENT USING INVALID BUT POSITIVE
C                  VALUES OF LINK THAT WILL PROBABLY EXTEND
C                  BEYOND THE PROGRAM LIMITS OF PRGOPT(*).
C                  UNRECOGNIZED VALUES OF KEY ARE IGNORED.  THE
C                  ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER
C                  OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING
C                  RESTRICTION.  TO PREVENT CYCLING IN THE
C                  PROCESSING OF THE OPTION ARRAY A COUNT OF THE
C                  NUMBER OF OPTIONS CHANGED IS MAINTAINED.
C                  WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR
C                  MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS.
C
C                  OPTIONS..
C
C                  KEY=1
C                         COMPUTE IN W(*,*) THE N BY N
C                  COVARIANCE MATRIX OF THE SOLUTION VARIABLES
C                  AS AN OUTPUT PARAMETER.  NOMINALLY THE
C                  COVARIANCE MATRIX WILL NOT BE COMPUTED.
C                  (THIS REQUIRES NO USER INPUT.)
C                  THE DATA SET FOR THIS OPTION IS A SINGLE VALUE.
C                  IT MUST BE NONZERO WHEN THE COVARIANCE MATRIX
C                  IS DESIRED.  IF IT IS ZERO, THE COVARIANCE
C                  MATRIX IS NOT COMPUTED.  WHEN THE COVARIANCE MATRIX
C                  IS COMPUTED, THE FIRST DIMENSIONING PARAMETER
C                  OF THE ARRAY W(*,*) MUST SATISFY MDW.GE.MAX0(M,N).
C
C                  KEY=2
C                         SCALE THE NONZERO COLUMNS OF THE
C                         ENTIRE DATA MATRIX.
C                  (E)
C                  (A)
C                  (G)
C
C                  TO HAVE LENGTH ONE.  THE DATA SET FOR THIS
C                  OPTION IS A SINGLE VALUE.  IT MUST BE
C                  NONZERO IF UNIT LENGTH COLUMN SCALING
C                  IS DESIRED.
C
C                  KEY=3
C                         SCALE COLUMNS OF THE ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  (G)
C
C                  WITH A USER-PROVIDED DIAGONAL MATRIX.
C                  THE DATA SET FOR THIS OPTION CONSISTS
C                  OF THE N DIAGONAL SCALING FACTORS, ONE FOR
C                  EACH MATRIX COLUMN.
C
C                  KEY=4
C                         CHANGE THE RANK DETERMINATION TOLERANCE FOR
C                  THE EQUALITY CONSTRAINT EQUATIONS FROM
C                  THE NOMINAL VALUE OF DSQRT(DRELPR).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN DRELPR, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY DRELPR IS THE
C                  LARGEST POSITIVE NUMBER SUCH THAT T=1.+DRELPR
C                  SATISFIES T.EQ.1.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST DRELPR.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  KEY=5
C                         CHANGE THE RANK DETERMINATION TOLERANCE FOR
C                  THE REDUCED LEAST SQUARES EQUATIONS FROM
C                  THE NOMINAL VALUE OF DSQRT(DRELPR).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN DRELPR, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST DRELPR.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  FOR EXAMPLE, SUPPOSE WE WANT TO CHANGE
C                  THE TOLERANCE FOR THE REDUCED LEAST SQUARES
C                  PROBLEM, COMPUTE THE COVARIANCE MATRIX OF
C                  THE SOLUTION PARAMETERS, AND PROVIDE
C                  COLUMN SCALING FOR THE DATA MATRIX.  FOR
C                  THESE OPTIONS THE DIMENSION OF PRGOPT(*)
C                  MUST BE AT LEAST N+9.  THE FORTRAN STATEMENTS
C                  DEFINING THESE OPTIONS WOULD BE AS FOLLOWS.
C
C                  PRGOPT(1)=4 (LINK TO ENTRY 4 IN PRGOPT(*))
C                  PRGOPT(2)=1 (COVARIANCE MATRIX KEY)
C                  PRGOPT(3)=1 (COVARIANCE MATRIX WANTED)
C
C                  PRGOPT(4)=7 (LINK TO ENTRY 7 IN PRGOPT(*))
C                  PRGOPT(5)=5 (LEAST SQUARES EQUAS. TOLERANCE KEY)
C                  PRGOPT(6)=... (NEW VALUE OF THE TOLERANCE)
C
C                  PRGOPT(7)=N+9 (LINK TO ENTRY N+9 IN PRGOPT(*))
C                  PRGOPT(8)=3 (USER-PROVIDED COLUMN SCALING KEY)
C
C                  CALL DCOPY(N,D,1,PRGOPT(9),1)  (COPY THE N
C                    SCALING FACTORS FROM THE USER ARRAY D(*)
C                    TO PRGOPT(9)-PRGOPT(N+8))
C
C                  PRGOPT(N+9)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  THE CONTENTS OF PRGOPT(*) ARE NOT MODIFIED
C                  BY THE SUBPROGRAM.
C                  THE OPTIONS FOR WNNLS( ) CAN ALSO BE INCLUDED
C                  IN THIS ARRAY.  THE VALUES OF KEY RECOGNIZED
C                  BY WNNLS( ) ARE 6, 7 AND 8.  THEIR FUNCTIONS
C                  ARE DOCUMENTED IN THE USAGE INSTRUCTIONS FOR
C                  SUBROUTINE WNNLS( ).  NORMALLY THESE OPTIONS
C                  DO NOT NEED TO BE MODIFIED WHEN USING LSEI( ).
C
C     IP(1),       THE AMOUNTS OF WORKING STORAGE ACTUALLY
C     IP(2)        ALLOCATED FOR THE WORKING ARRAYS WS(*) AND
C                  IP(*), RESPECTIVELY.  THESE QUANTITIES ARE
C                  COMPARED WITH THE ACTUAL AMOUNTS OF STORAGE
C                  NEEDED BY LSEI( ).  INSUFFICIENT STORAGE
C                  ALLOCATED FOR EITHER WS(*) OR IP(*) IS AN
C                  ERROR.  THIS FEATURE WAS INCLUDED IN LSEI( )
C                  BECAUSE MISCALCULATING THE STORAGE FORMULAS
C                  FOR WS(*) AND IP(*) MIGHT VERY WELL LEAD TO
C                  SUBTLE AND HARD-TO-FIND EXECUTION ERRORS.
C
C                  THE LENGTH OF WS(*) MUST BE AT LEAST
C
C                  LW = 2*(ME+N)+K+(MG+2)*(N+7)
C
C
C                  WHERE K = MAX(MA+MG,N)
C                  THIS TEST WILL NOT BE MADE IF IP(1).LE.0.
C
C                  THE LENGTH OF IP(*) MUST BE AT LEAST
C
C                  LIP = MG+2*N+2
C                  THIS TEST WILL NOT BE MADE IF IP(2).LE.0.
C
C     OUTPUT..
C
C     X(*),RNORME,  THE ARRAY X(*) CONTAINS THE SOLUTION PARAMETERS
C     RNORML        IF THE INTEGER OUTPUT FLAG MODE = 0 OR 1.
C                   THE DEFINITION OF MODE IS GIVEN DIRECTLY BELOW.
C                   WHEN MODE = 0 OR 1, RNORME AND RNORML
C                   RESPECTIVELY CONTAIN THE RESIDUAL VECTOR
C                   EUCLIDEAN LENGTHS OF F - EX AND B - AX.  WHEN
C                   MODE=1 THE EQUALITY CONSTRAINT EQUATIONS EX=F
C                   ARE CONTRADICTORY, SO RNORME.NE.0. THE RESIDUAL
C                   VECTOR F-EX HAS MINIMAL EUCLIDEAN LENGTH. FOR
C                   MODE.GE.2, NONE OF THESE PARAMETERS ARE
C                   DEFINED.
C
C     MODE          INTEGER FLAG THAT INDICATES THE SUBPROGRAM
C                   STATUS AFTER COMPLETION.  IF MODE.GE.2, NO
C                   SOLUTION HAS BEEN COMPUTED.
C
C                   MODE =
C
C                   0  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
C                      ARE COMPATIBLE AND HAVE BEEN SATISFIED.
C
C                   1  EQUALITY CONSTRAINTS ARE CONTRADICTORY.
C                      A GENERALIZED INVERSE SOLUTION OF EX=F WAS USED
C                      TO MINIMIZE THE RESIDUAL VECTOR LENGTH F-EX.
C                      IN THIS SENSE, THE SOLUTION IS STILL MEANINGFUL.
C
C                   2  INEQUALITY CONSTRAINTS ARE CONTRADICTORY.
C
C                   3  BOTH EQUALITY AND INEQUALITY CONSTRAINTS
C                      ARE CONTRADICTORY.
C
C                   THE FOLLOWING INTERPRETATION OF
C                   MODE=1,2 OR 3 MUST BE MADE.  THE
C                   SETS CONSISTING OF ALL SOLUTIONS
C                   OF THE EQUALITY CONSTRAINTS EX=F
C                   AND ALL VECTORS SATISFYING GX.GE.H
C                   HAVE NO POINTS ON COMMON.  (IN
C                   PARTICULAR THIS DOES NOT SAY THAT
C                   EACH INDIVIDUAL SET HAS NO POINTS
C                   AT ALL, ALTHOUGH THIS COULD BE THE
C                   CASE.)
C
C                   4  USAGE ERROR OCCURRED.  THE VALUE
C                      OF MDW IS .LT. ME+MA+MG, MDW IS
C                      .LT. N AND A COVARIANCE MATRIX IS
C                      REQUESTED, THE OPTION VECTOR
C                      PRGOPT(*) IS NOT PROPERLY DEFINED,
C                      OR THE LENGTHS OF THE WORKING ARRAYS
C                      WS(*) AND IP(*), WHEN SPECIFIED IN
C                      IP(1) AND IP(2) RESPECTIVELY, ARE NOT
C                      LONG ENOUGH.
C
C     W(*,*)        THE ARRAY W(*,*) CONTAINS THE N BY N SYMMETRIC
C                   COVARIANCE MATRIX OF THE SOLUTION PARAMETERS,
C                   PROVIDED THIS WAS REQUESTED ON INPUT WITH
C                   THE OPTION VECTOR PRGOPT(*) AND THE OUTPUT
C                   FLAG IS RETURNED WITH MODE = 0 OR 1.
C
C     IP(*)         THE INTEGER WORKING ARRAY HAS THREE ENTRIES
C                   THAT PROVIDE RANK AND WORKING ARRAY LENGTH
C                   INFORMATION AFTER COMPLETION.
C
C                      IP(1) = RANK OF EQUALITY CONSTRAINT
C                              MATRIX.  DEFINE THIS QUANTITY
C                              AS KRANKE.
C
C                      IP(2) = RANK OF REDUCED LEAST SQUARES
C                              PROBLEM.
C
C                      IP(3) = THE AMOUNT OF STORAGE IN THE
C                              WORKING ARRAY WS(*) THAT WAS
C                              ACTUALLY USED BY THE SUBPROGRAM.
C                              THE FORMULA GIVEN ABOVE FOR THE LENGTH
C                              OF WS(*) IS A NECESSARY OVERESTIMATE.
C     USER DESIGNATED
C     WORKING ARRAYS..
C
C     WS(*),IP(*)              THESE ARE RESP. TYPE FLOATING POINT
C                              AND TYPE INTEGER WORKING ARRAYS.
C                              THEIR REQUIRED MINIMAL LENGTHS ARE
C                              GIVEN ABOVE.
C
C
C     SUBROUTINES CALLED
C
C     LSI           PART OF THIS PACKAGE.  SOLVES A
C                   CONSTRAINED LEAST SQUARES PROBLEM WITH
C                   INEQUALITY CONSTRAINTS.
C
C++
C     DDOT,DSCAL,   SUBROUTINES FROM THE BLAS PACKAGE.
C     DAXPY,DASUM,  SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308.
C     DCOPY,DNRM2,
C     DSWAP
C
C     H12           SUBROUTINE TO CONSTRUCT AND APPLY A
C                   HOUSEHOLDER TRANSFORMATION.
C
C     XERROR        FROM SLATEC ERROR PROCESSING PACKAGE.
C                   THIS IS DOCUMENTED IN SANDIA TECH. REPT.,
C                   SAND78-1189.
C
C     SUBROUTINE LSEI(W,MDW,ME,MA,MG,N,PRGOPT,X,RNORME,RNORML,MODE,WS,
C    1 IP)
C
C     REVISED OCT. 1, 1981.
C
      DOUBLE PRECISION W(MDW,1), PRGOPT(1), X(1), WS(1), RNORME, RNORML
      INTEGER IP(3)
      DOUBLE PRECISION DUMMY, ENORM, DRELPR, FNORM, GAM, HALF, ONE, RB,
     *RN, RNMAX, SIZE, SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME, ZERO
      DOUBLE PRECISION DASUM, DDOT, DNRM2, DSQRT, DABS, DMAX1
      LOGICAL COV
      DATA ZERO /0.D0/, DRELPR /0.D0/, HALF /0.5E0/, ONE /1.D0/
C
C     CHECK THAT ENOUGH STORAGE WAS ALLOCATED IN WS(*) AND IP(*).
      IF (.NOT.(IP(1).GT.0)) GO TO 20
      LCHK = 2*(ME+N) + MAX0(MA+MG,N) + (MG+2)*(N+7)
      IF (.NOT.(IP(1).LT.LCHK)) GO TO 10
      MODE = 4
      NERR = 2
      IOPT = 1
      CALL XERRWV(67HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR WS(*),
     *NEED LW=I1 BELOW, 67, NERR, IOPT, 1, LCHK, 0,
     * 0, DUMMY, DUMMY)
      RETURN
   10 CONTINUE
   20 IF (.NOT.(IP(2).GT.0)) GO TO 40
      LCHK = MG + 2*N + 2
      IF (.NOT.(IP(2).LT.LCHK)) GO TO 30
      MODE = 4
      NERR = 2
      IOPT = 1
      CALL XERRWV(68HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR IP(*),
     *NEED LIP=I1 BELOW, 68, NERR, IOPT, 1, LCHK, 0,
     * 0, DUMMY, DUMMY)
      RETURN
   30 CONTINUE
C
C     COMPUTE MACHINE PRECISION=DRELPR ONLY WHEN NECESSARY.
   40 IF (.NOT.(DRELPR.EQ.ZERO)) GO TO 70
      DRELPR = ONE
   50 IF (ONE+DRELPR.EQ.ONE) GO TO 60
      DRELPR = DRELPR*HALF
      GO TO 50
   60 DRELPR = DRELPR + DRELPR
C
C     COMPUTE NUMBER OF POSSIBLE RIGHT MULTIPLYING HOUSEHOLDER
C     TRANSFORMATIONS.
   70 M = ME + MA + MG
      MODE = 0
      IF (N.LE.0 .OR. M.LE.0) RETURN
      IF (.NOT.(MDW.LT.M)) GO TO 80
      NERR = 1
      IOPT = 1
      CALL XERROR(36HLSEI( ), MDW.LT.ME+MA+MG IS AN ERROR, 36, NERR,
     * IOPT)
      MODE = 4
      RETURN
   80 NP1 = N + 1
      KRANKE = MIN0(ME,N)
      N1 = 2*KRANKE + 1
      N2 = N1 + N
C
C     PROCESS-OPTION-VECTOR
      ASSIGN 90 TO IGO990
      GO TO 480
   90 IF (.NOT.(COV .AND. MDW.LT.N)) GO TO 100
      NERR = 2
      IOPT = 1
      CALL XERROR(
     * 54HLSEI( ), MDW.LT.N, WHEN COV MATRIX NEEDED, IS AN ERROR, 54,
     * NERR, IOPT)
      MODE = 4
      RETURN
  100 L = KRANKE
C
C     COMPUTE NORM OF EQUALITY CONSTRAINT MATRIX AND RT SIDE.
      ENORM = ZERO
      DO 110 J=1,N
        ENORM = DMAX1(ENORM,DASUM(ME,W(1,J),1))
  110 CONTINUE
      FNORM = DASUM(ME,W(1,NP1),1)
      IF (.NOT.(L.GT.0)) GO TO 190
      SNMAX = ZERO
      RNMAX = ZERO
      DO 180 I=1,L
C
C     COMPUTE MAXIMUM RATIO OF VECTOR LENGTHS. PARTITION
C     IS AT COL. I.
        DO 150 K=I,ME
          SN = DDOT(N-I+1,W(K,I),MDW,W(K,I),MDW)
          RN = DDOT(I-1,W(K,1),MDW,W(K,1),MDW)
          IF (.NOT.(RN.EQ.ZERO .AND. SN.GT.SNMAX)) GO TO 120
          SNMAX = SN
          IMAX = K
          GO TO 140
  120     IF (.NOT.(K.EQ.I .OR. (SN*RNMAX.GT.RN*SNMAX))) GO TO 130
          SNMAX = SN
          RNMAX = RN
          IMAX = K
  130     CONTINUE
  140     CONTINUE
  150   CONTINUE
C
C     INTERCHANGE ROWS IF NECESSARY.
        IF (I.NE.IMAX) CALL DSWAP(NP1, W(I,1), MDW, W(IMAX,1), MDW)
        IF (.NOT.(SNMAX.GT.TAU**2*RNMAX)) GO TO 160
C
C     ELIMINATE ELEMS I+1,...,N IN ROW I.
        CALL H12(1, I, I+1, N, W(I,1), MDW, WS(I), W(I+1,1), MDW, 1,
     *   M-I)
        GO TO 170
  160   KRANKE = I - 1
        GO TO 200
  170   CONTINUE
  180 CONTINUE
  190 CONTINUE
  200 CONTINUE
C
C     SAVE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL DCOPY(KRANKE, W, MDW+1, WS(KRANKE+1), 1)
C
C     USE HOUSEHOLDER TRANS FROM LEFT TO ACHIEVE KRANKE BY KRANKE UPPER
C     TRIANGULAR FORM.
      IF (.NOT.(KRANKE.GT.0 .AND. KRANKE.LT.ME)) GO TO 220
      DO 210 KK=1,KRANKE
        K = KRANKE + 1 - KK
C
C     APPLY TRANFORMATION TO MATRIX COLS. 1,...,K-1.
        CALL H12(1, K, KRANKE+1, ME, W(1,K), 1, UP, W, 1, MDW, K-1)
C
C     APPLY TO RT SIDE VECTOR.
        CALL H12(2, K, KRANKE+1, ME, W(1,K), 1, UP, W(1,NP1), 1, 1, 1)
  210 CONTINUE
  220 IF (.NOT.(KRANKE.GT.0)) GO TO 240
C
C     SOLVE FOR VARIABLES 1,...,KRANKE IN NEW COORDINATES.
      CALL DCOPY(KRANKE, W(1,NP1), 1, X, 1)
      DO 230 I=1,KRANKE
        X(I) = (X(I)-DDOT(I-1,W(I,1),MDW,X,1))/W(I,I)
  230 CONTINUE
C
C     COMPUTE RESIDUALS FOR REDUCED PROBLEM.
  240 MEP1 = ME + 1
      RNORML = ZERO
      IF (.NOT.(ME.LT.M)) GO TO 270
      DO 260 I=MEP1,M
        W(I,NP1) = W(I,NP1) - DDOT(KRANKE,W(I,1),MDW,X,1)
        SN = DDOT(KRANKE,W(I,1),MDW,W(I,1),MDW)
        RN = DDOT(N-KRANKE,W(I,KRANKE+1),MDW,W(I,KRANKE+1),MDW)
        IF (.NOT.(RN.LE.TAU**2*SN .AND. KRANKE.LT.N)) GO TO 250
        W(I,KRANKE+1) = ZERO
        CALL DCOPY(N-KRANKE, W(I,KRANKE+1), 0, W(I,KRANKE+1), MDW)
  250   CONTINUE
  260 CONTINUE
C
C     COMPUTE EQUAL. CONSTRAINT EQUAS. RESIDUAL LENGTH.
  270 RNORME = DNRM2(ME-KRANKE,W(KRANKE+1,NP1),1)
C
C     MOVE REDUCED PROBLEM DATA UPWARD IF KRANKE.LT.ME.
      IF (.NOT.(KRANKE.LT.ME)) GO TO 290
      DO 280 J=1,NP1
        CALL DCOPY(M-ME, W(ME+1,J), 1, W(KRANKE+1,J), 1)
  280 CONTINUE
C
C     COMPUTE SOLN OF REDUCED PROBLEM.
  290 CALL LSI(W(KRANKE+1,KRANKE+1), MDW, MA, MG, N-KRANKE, PRGOPT,
     * X(KRANKE+1), RNORML, MODE, WS(N2), IP(2))
      IF (.NOT.(ME.GT.0)) GO TO 330
C
C     TEST FOR CONSISTENCY OF EQUALITY CONSTRAINTS.
      MDEQC = 0
      XNRME = DASUM(KRANKE,W(1,NP1),1)
      IF (RNORME.GT.TAU*(ENORM*XNRME+FNORM)) MDEQC = 1
      MODE = MODE + MDEQC
C
C     CHECK IF SOLN TO EQUAL. CONSTRAINTS SATISFIES INEQUAL.
C     CONSTRAINTS WHEN THERE ARE NO DEGREES OF FREEDOM LEFT.
      IF (.NOT.(KRANKE.EQ.N .AND. MG.GT.0)) GO TO 320
      XNORM = DASUM(N,X,1)
      MAPKE1 = MA + KRANKE + 1
      MEND = MA + KRANKE + MG
      DO 310 I=MAPKE1,MEND
        SIZE = DASUM(N,W(I,1),MDW)*XNORM + DABS(W(I,NP1))
        IF (.NOT.(W(I,NP1).GT.TAU*SIZE)) GO TO 300
        MODE = MODE + 2
        GO TO 450
  300   CONTINUE
  310 CONTINUE
  320 CONTINUE
  330 IF (.NOT.(KRANKE.GT.0)) GO TO 420
C
C     REPLACE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL DCOPY(KRANKE, WS(KRANKE+1), 1, W, MDW+1)
C
C     REAPPLY TRANS TO PUT SOLN IN ORIGINAL COORDINATES.
      DO 340 II=1,KRANKE
        I = KRANKE + 1 - II
        CALL H12(2, I, I+1, N, W(I,1), MDW, WS(I), X, 1, 1, 1)
  340 CONTINUE
C
C     COMPUTE COV MATRIX OF EQUAL. CONSTRAINED PROBLEM.
      IF (.NOT.(COV)) GO TO 410
      DO 400 JJ=1,KRANKE
        J = KRANKE + 1 - JJ
        IF (.NOT.(J.LT.N)) GO TO 390
        RB = WS(J)*W(J,J)
        IF (RB.NE.ZERO) RB = ONE/RB
        JP1 = J + 1
        DO 350 I=JP1,N
          W(I,J) = DDOT(N-J,W(I,JP1),MDW,W(J,JP1),MDW)*RB
  350   CONTINUE
        GAM = DDOT(N-J,W(JP1,J),1,W(J,JP1),MDW)*RB
        GAM = HALF*GAM
        CALL DAXPY(N-J, GAM, W(J,JP1), MDW, W(JP1,J), 1)
        DO 370 I=JP1,N
          DO 360 K=I,N
            W(I,K) = W(I,K) + W(J,I)*W(K,J) + W(I,J)*W(J,K)
            W(K,I) = W(I,K)
  360     CONTINUE
  370   CONTINUE
        UJ = WS(J)
        VJ = GAM*UJ
        W(J,J) = UJ*VJ + UJ*VJ
        DO 380 I=JP1,N
          W(J,I) = UJ*W(I,J) + VJ*W(J,I)
  380   CONTINUE
        CALL DCOPY(N-J, W(J,JP1), MDW, W(JP1,J), 1)
  390   CONTINUE
  400 CONTINUE
  410 CONTINUE
C
C     APPLY THE SCALING TO THE COVARIANCE MATRIX.
  420 IF (.NOT.(COV)) GO TO 440
      DO 430 I=1,N
        L = N1 + I
        CALL DSCAL(N, WS(L-1), W(I,1), MDW)
        CALL DSCAL(N, WS(L-1), W(1,I), 1)
  430 CONTINUE
  440 CONTINUE
  450 CONTINUE
C
C     RESCALE SOLN. VECTOR.
      IF (.NOT.(MODE.LE.1)) GO TO 470
      DO 460 J=1,N
        L = N1 + J
        X(J) = X(J)*WS(L-1)
  460 CONTINUE
  470 IP(1) = KRANKE
      IP(3) = IP(3) + 2*KRANKE + N
      RETURN
  480 CONTINUE
C     TO PROCESS-OPTION-VECTOR
C
C     THE NOMINAL TOLERANCE USED IN THE CODE
C     FOR THE EQUALITY CONSTRAINT EQUATIONS.
      TAU = DSQRT(DRELPR)
C
C     THE NOMINAL COLUMN SCALING USED IN THE CODE IS
C     THE IDENTITY SCALING.
      WS(N1) = ONE
      CALL DCOPY(N, WS(N1), 0, WS(N1), 1)
C
C     NO COVARIANCE MATRIX IS NOMINALLY COMPUTED.
      COV = .FALSE.
C
C     DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE.
      NOPT = 1000
      NTIMES = 0
C
C     DEFINE BOUND FOR POSITIVE VALUES OF LINK.
      NLINK = 100000
      LAST = 1
      LINK = PRGOPT(1)
      IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 490
      NERR = 3
      IOPT = 1
      CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
      MODE = 4
      RETURN
  490 IF (.NOT.(LINK.GT.1)) GO TO 540
      NTIMES = NTIMES + 1
      IF (.NOT.(NTIMES.GT.NOPT)) GO TO 500
      NERR = 3
      IOPT = 1
      CALL XERROR(
     * 52HLSEI( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 52,
     * NERR, IOPT)
      MODE = 4
      RETURN
  500 KEY = PRGOPT(LAST+1)
      IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO
      IF (.NOT.(KEY.EQ.2 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 520
      DO 510 J=1,N
        T = DNRM2(M,W(1,J),1)
        IF (T.NE.ZERO) T = ONE/T
        L = N1 + J
        WS(L-1) = T
  510 CONTINUE
  520 IF (KEY.EQ.3) CALL DCOPY(N, PRGOPT(LAST+2), 1, WS(N1), 1)
      IF (KEY.EQ.4) TAU = DMAX1(DRELPR,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 530
      NERR = 3
      IOPT = 1
      CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
      MODE = 4
      RETURN
  530 LAST = LINK
      LINK = NEXT
      GO TO 490
  540 DO 550 J=1,N
        L = N1 + J
        CALL DSCAL(M, WS(L-1), W(1,J), 1)
  550 CONTINUE
      GO TO 560
  560 GO TO IGO990, (90)
      END
